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ABSTRACT 

Using a linear stability analysis and two and three-dimensional nonlinear simulations, we study 
the physics of buoyancy instabilities in a combined thermal and relativistic (cosmic ray) plasma, 
motivated by the application to clusters of galaxies. We argue that cosmic ray diffusion is likely to 
be slow compared to the buoyancy time on large length scales, so that cosmic rays are effectively 
adiabatic. If the cosmic ray pressure p cr is > 25% of the thermal pressure, and the cosmic ray 
entropy (p C r/p 4 ^ 3 ', P is the thermal plasma density) decreases outwards, cosmic rays drive an adiabatic 
convective instability analogous to Schwarzschild convection in stars. Global simulations of galaxy 
cluster cores show that this instability saturates by reducing the cosmic ray entropy gradient and 
driving efficient convection and turbulent mixing. At larger radii in cluster cores, the thermal plasma 
is unstable to the heat flux-driven buoyancy instability (HBI), a convective instability generated 
by anisotropic thermal conduction and a background conductive heat flux. The HBI saturates by 
rearranging the magnetic field lines to become largely perpendicular to the local gravitational field; 
the resulting turbulence also primarily mixes plasma in the perpendicular plane. Cosmic-ray driven 
convection and the HBI may contribute to redistributing metals produced by Type la supernovae in 
clusters. Our calculations demonstrate that adiabatic simulations of galaxy clusters can artificially 
suppress the mixing of thermal and relativistic plasma; anisotropic thermal conduction allows more 
efficient mixing, which may contribute to cosmic rays being distributed throughout the cluster volume. 
Subject headings: convection — cooling flows — galaxies: active — galaxies: clusters: general - 
magnetic fields 



1. INTRODUCTION 

Microscopic transport of heat and momentum in di- 
lute plasmas, like those in clusters of galaxies, is primar- 
ily along magnetic field lines (Braginskii 1965). This 
anisotropic transport dramatically affects the convective 
stability of the plasma; convective stability is no longer 
determined by the entropy gradient (Schwarzschild 
1958). Instead, a plasma is unstable to buoyant mo- 
tions irrespective of the background entropy and temper- 
ature gradients (Balbus 2000; Quataert 2008). Cosmic 
rays diffusing along magnetic field lines also affect the 
convective stability of the plasma (Chandran & Dennis 
2006; Dennis & Chandran 2008). Nonlinear simulations 
show that these instabilities driven by anisotropic ther- 
mal and cosmic ray transport can change the magnetic 
field configuration, and the background temperature and 
density profiles in the plasma, but they do not drive ef- 
ficient convection (e.g., Parrish & Stone 2007; Parrish 
& Quataert 2008; Parrish, Stone, & Lemaster 2008; 
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Sharma, Quataert, & Stone 2008). The instabilities sat- 
urate largely by rearranging the magnetic field configu- 
ration, thereby slowing down the instability and reaching 
a state of marginal stability to linear perturbations. By 
contrast, hydrodynamic convection in stellar interiors re- 
distributes energy efficiently to make the plasma nearly 
adiabatic and thus marginally stable to convection. 

One of the key astrophysical motivations for studying 
the transport properties of dilute plasmas in the pres- 
ence of cosmic rays is to understand the dynamical and 
thermal structure of clusters of galaxies. The radiative 
cooling time (< 1 Gyr) is much less than the Hubble 
time (tn ~ 13.7 Gyr) in cluster cores. Thus, it was ex- 
pected that the intracluster medium (ICM) would cool 
rapidly, resulting in large rates (> 100M Q yr" 1 ) of mass 
cooling to form cold gas and stars (e.g., Fabian 1994). 
However, X-ray observations have failed to detect copi- 
ous emission from the expected cold plasma component 
in cluster cores (e.g., Peterson et al. 2003). The lack 
of cooling flows implies that cooling is balanced by some 
source of heating, e.g., heating by thermal conduction 
from large radii (Bertschinger & Meiksin 1986; Zakamska 
& Narayan 2003), heating by jets and bubbles blown by 
a central AGN (Binney & Tabor 1995; Ciotti & Ostriker 
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2001), or heating by cosmic rays (e.g., Rosner & Tucker 
1989; Loewenstein, Zweibel, & Begelman 1991; Chan- 
dran & Rasera 2007). Although thermal conduction 
may operate at large radii, it appears that the plasma at 
small radii must be heated by a feedback process which 
efficiently self-regulates. This is required to avoid the fine 
tuning of thermal conductivity required in models that 
include only conduction (e.g., Guo & Oh 2008; Conroy 
& Ostriker' 2008). 

Cosmic rays from a central AGN have been invoked 
to prevent catastrophic cooling of the plasma in cluster 
cores, either directly via Alfven waves driven by cosmic 
rays heating the plasma (e.g., Loewenstein, Zweibel, & 
Begelman 1991; Guo & Oh 2008), or indirectly via con- 
vection driven by cosmic rays, the dissipation of which 
heats the plasma (Chandran & Rasera 2007; Rebusco 
et al. 2006). Although there is ample evidence for the 
presence of cosmic rays in radio emitting bubbles in clus- 
ters (e.g., Birzan et al. 2004), it is unclear how/whether 
cosmic rays can be spread throughout the cluster vol- 
ume at a sufficient level for these heating mechanisms to 
work. Simple hydrodynamic jets do not couple their en- 
ergy to most of the ICM and instead simply drill through 
it, without heating and without transporting cosmic rays 
(and metals) throughout the ICM (Vernaleo & Reynolds 
2006). 

In this paper we show that cosmic rays, which are 
likely to be centrally concentrated in clusters, can drive 
efficient convection and mixing if the cosmic ray pres- 
sure is not negligible compared to the plasma pressure 
{Pcr/p ^ 0.25); this is likely the case in and around radio 
bubbles. We argue that on the large scales that likely 
dominate the turbulent dynamics in the ICM, cosmic 
rays are effectively adiabatic rather than diffusive (i.e., 
the cosmic ray diffusion time is longer than the buoy- 
ancy time). As a result, the cosmic rays can drive a 
Schwarzschild-like adiabatic convective instability. We 
present a linear analysis demonstrating that, while mag- 
netic reorientation can shut off diffusive ( "isobaric" ) cos- 
mic ray instabilities, it cannot shut off the adiabatic 
buoyancy instability driven by a negative cosmic ray en- 
tropy gradient. We then present two- fluid (plasma and 
cosmic rays) numerical simulations with thermal conduc- 
tion and cosmic ray diffusion along magnetic field lines. 

We do not include plasma cooling in this paper, nor do 
we include the heating of the thermal plasma that arises 
from the excitation of short wavelength Alfven waves by 
streaming cosmic rays. These are both significant omis- 
sions and preclude our results from being an accurate 
representation of the plasma in cluster cores. However, 
the main focus of this paper is not to solve the cooling 
flow problem per se, but rather to isolate and under- 
stand the transport and turbulence properties of cluster 
plasmas with realistic physics (e.g., anisotropic conduc- 
tion, convection, and cosmic rays). By neglecting cool- 
ing, our calculations implicitly assume that some unspec- 
ified source of heating is preventing the rapid cooling 
of the ICM. A study of cluster cores with cooling and 
anisotropic conduction will be presented in a separate 
paper. 

The remainder of this paper is organized as follows. 
In §2 we present the basic equations used in our anal- 
ysis and derive the dispersion relation for linear buoy- 
ancy waves in the presence of thermal plasma and cosmic 



rays. In §3 we describe our numerical simulations and 
show that a steep cosmic-ray entropy gradient drives con- 
vective motions at small radii in cluster cores, and that 
anisotropic thermal conduction drives convection at in- 
termediate radii that rearranges the magnetic field struc- 
ture in clusters. Readers interested in just the numeri- 
cal results can skip the linear stability calculation in §2. 
In §4 we summarize and discuss the implications of our 
work. 

2. BASIC EQUATIONS 

To describe cosmic rays and thermal plasma in the 
ICM, we use the two-fluid model of Drury & Volk (1981), 
modified to include anisotropic transport, gravitational 
acceleration g = —gr (g = d$/dr, where $ is the gravi- 
tational potential), and a cosmic ray source term to build 
up cosmic ray pressure. The equations of this model are 



dp 
~dt 



= -pV • v, 



(1) 



dv (V x B) x B 
P-JI = - 7Z V( P + p CI )- pgr, (2) 



dt 



47T 



dB 

~dt 



= V x (v x B), 



§-^§ = -(7-DV.Q, 
dt p dt 



(3) 



(4) 



and 



dp cr 7crPcr dp , . . . 

where d/dt = d/dt + v ■ V is the Lagrangian time deriva- 
tive, 

Q = —K\\b(b • VT) (6) 
is the heat flux, T is the plasma temperature, 



r= -Dub(b- Vp cr ) 



(7) 



is the diffusive flux of cosmic-ray energy (multiplied by 
[7cr — 1] ) ! Qc is the cosmic ray energy source term, p is the 
mass density, v is the common bulk-flow velocity of the 
thermal plasma and cosmic rays, B is the magnetic field, 
b = B I B, p and p cr are the thermal-plasma and cosmic- 
ray pressures, k\\ is the parallel thermal conductivity, D|| 
is the diffusion coefficient for cosmic-ray transport along 
the magnetic field, and 7=5/3 and 7 cr =4/3 are the adi- 
abatic indices of the thermal plasma and cosmic rays, 
respectively. 

As mentioned in §1, we do not include radiative cool- 
ing in the energy equation (eq. [4]). In addition, we do 
not include the effects of cosmic ray streaming relative to 
the thermal plasma: in particular we neglect Alfven wave 
heating of the thermal plasma and Alfven wave stream- 
ing in the cosmic ray energy equation (e.g., Loewenstein, 
Zweibel, & Begelman 1991). This physics will be in- 
cluded in the future together with plasma cooling. For 
the present paper we focus on the physics of convective 
instabilities in clusters using a simplified but physically 
reasonable model. 



2.1. Linear Stability Analysis 

We take all quantities to be the sum of an equilibrium 
value plus a small-amplitude fluctuation: B = Bq + B\, 
etc. We take the equilibrium velocity Vq to vanish, set 
Bq = Bq x x+Bo z z (z = r for the cluster, and x is chosen 
such that local magnetic field lies in the x — z plane), 
and take To, p C rO, and po to be functions of z alone. We 
employ a local analysis in which all fluctuating quantities 

vary as e l ^' x ~ lujt with kH 3> 1, where H is the scale 
on which the equilibrium quantities vary. We consider 
the limit (3 = 8npo/B 2 ;§> 1 and work in the Boussinesq 
approximation, lo <C kc s (c s is the sound speed). We also 
do not include the perturbed cosmic ray source term (eq. 
[5]) in our linear analysis since its form is uncertain. 

In terms of the plasma displacement, £ = — , the per- 
turbed magnetic field [B\ = z[fc- J3 ]£) can be combined 
with equation (2), to give 



Polo £ = 



(k-B Q ) 2 Z 
4ir 



+ ifclli + pigz, 



(8) 



where U\ = pi + p CI i + Bo ■ Bi/4it is the total-pressure 
perturbation. Dotting equation (8) with fc and using the 
near-incompressibility condition (fe • £ ~ £/H), we find 
to leading order in (kH)^ 1 that 



111 = ipigk z /k 2 



and 



pig ( „ kk z 



(w -k ]l v A )Z = — [z- — 



(9) 
(10) 



where fc|| = fc • 6 , and v\ = B^/A-upo is the square of the 
Alfven speed. 

In terms of the perturbation to the magnetic-field unit 
vector, 6i(= Bi/B - b b ■ Bi/B ) = ifc||[£ - b (b ■ 
£)], the perturbed heat flux can be written as Q l — 

-k\\ , [6i (bo ■ VT ) + b (&i ■ VT ) + 6 (So • VTi )] , where we 
have dropped the term involving the perturbed thermal 
conductivity since it is smaller than the other terms by a 
factor of <~ (kH) -1 . Taking the dot product of the heat 
flux with fc we find that 



-ik ■ Q 1 = K|| >0 fcf(2£||& 2 



6)^-*|«ii,ori, (ii) 



where Z\\ = Z'^o and b z is the z component of bo- Using 
equations (4) and (11), and using pi/po = Pi/ Po + Ti/Tq, 
we find that 



Pi- 



+ 



jLop Q N 2 £ z 
(lo + iv)g 
(jlo + iv)popi 
(uj + iv)po 



ivp (2£\\b z - Zz) 



uj + iv 



dhiTo 
dz 



(12) 



where v = (7 — l)fcj?K|| .o^o/Po is the ra te at which ther- 
mal conductivity smoothes out temperature fluctuations 
along the magnetic field, and 

N 2 = {gh)j-Hpo/pl) 
is the square of the Brunt- Vaisala frequency. 



In the same way that we obtained equations (11) 
and (12), for the cosmic rays we find that 



Perl ; 



j cr uip cr0 M 2 £ z 
(uj + irf)g 

Tcr^PcrOPl 



uj + if] 



dp, 



crO 



dz 



(lo + ir])p ' 

where r\ = k 2 D\\ is the rate at which diffusion smoothes 
out variations in p CI along the magnetic field, and 

M 2 = (3/7cr)^ln(pcroMn 

is the square of the Brunt- Vaisala frequency associated 
with the cosmic ray pressure. 

Adding equations (12) and (13), making use of equa- 
tion (9), and noting that the right-hand side of equa- 
tion (9) is much smaller than the individual terms on 
the right-hand sides of equations (12) and (13), we find 
that 



Po 
where 



(14) 



N =^ujN 2 /(uj + w)+j cr aujM 2 /(uj + iri), 



a = Pcro/PO 

u\ = g 
and 



lo + if 



d\nT 

dz 



+ 



iar) \ d\np a0 



uj + if] 



dz 



S = ('JLO + iv)/(uj + iv) + aj cr uj/ (uj + it]). 

Upon substituting equation (14) into equation (10), we 
obtain an equation for the plasma displacement alone, 



2jfe||£||ff 



+ 6 N + lo((£ z - 2^b z 



0. 



(15) 



We set £ = CiCi + £,2&2, where ei = fc x z/\k x z\ and 
e 2 = fc x e\. After taking the dot product of equa- 
tion (15) with ei and e 2 , we obtain two equations which 
can be written in matrix form as 



P-71 2 

I A 



An 





A22 



(16) 



Setting the determinant of the matrix on the left-hand 
side of equation (16) equal to zero, we obtain the dis- 
persion relation (lo 2 — k 2 v 2 A )A 2 2 — 0, where 2 A22 = 



i'7'-\ - 5 x sin 2 (^N 2 + JujfJ, 9 is the angle be- 



k 2 v 2 



tween fc and z, 

J=l-2b 2 z + 2b x b z k x k z /k 2 ± , 

2 We have dropped a term on the right-hand side of expression 
for A22 equal to — 2i<5~ 1 sin 2 Ok^gft^ 1 (bxkxkzk^ 2 — b z ), which is 
small for (i ;P kH; in the opposite limit ui = knVA + *7 with 



j/u> ~ (kH)- 
kH/3- 1 : 1. 



Notice that k^v 2 A 



kug/3- 1 : N :: (kH) 2 ft- 1 
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is a factor that depends on the magnetic field geome- 
try and wavenumber, k\ — k 2 + k 2 7 and b x is the x- 

component of b . 
While one solution to equation (16), ui 2 = 

scribes waves unaffected by buoyancy, A22 
sponds to the modes modified by buoyancy, 



1.2 2 

«,,, u A , 



de- 



corre- 



ct 



- fcft£ - (T 1 sin 2 (jV 2 + Jw 2 ) =0. (17) 

In the limit that p a /p —* 0, equation (17) reduces to 
equation (13) of Quataert (2008). For nonzero a {p C r/p), 
we consider two limiting cases of equation (17). The 
first is a highly diffusive limit for cosmic rays, in which 
{77, v} ^> {oj,N, Ldi} ;§> k\\VA- In this case, equation (17) 
reduces to 



oj =gj 



din 7} 







dz 



1 dp cA 
dz 



Po 



(18) 



This is a generalization of the HBI (heat-flux driven 
buoyancy instability; Quataert 2008) and MTI (mag- 
netothermal instability; Balbus 2000) in the presence 
of diffusive cosmic rays (e.g., Chandran & Dennis 2006; 
Dennis & Chandran ' 2008). 

The second limit we consider is that of adiabatic cosmic 
rays, with v 3> {w, N, u>i} 3> {77, AimWa}- In this case, 
the cosmic rays diffuse slowly compared to the buoyancy 
time and behave nearly adiabatically; equation (17) then 
reduces to 



9 q sin 9 



1 + «7c 



PcrO d ( P„o \ 
po dz \Po" J 



J^lnTo 

dz 



■ (19) 



This implies that if the cosmic ray pressure is significant 
and if cosmic ray entropy (p C rO / Po" ) is sufficiently peaked 
and decreasing outward, the plasma will become unsta- 
ble to a Schwarzschild type buoyancy instability. The 
entropy of cluster plasma is stably stratified according 
to the Schwarzschild criterion. The thermal plasma re- 
sponse at small scales is nonetheless governed by the tem- 
perature gradient and not the entropy gradient. For typ- 
ical cluster parameters, although the global conduction 
timescale may be longer than the buoyancy timescale 
(see Fig. 1), anisotropic conduction determines the local 
buoyant response (i.e., v > N for kr 3> 1) via HBI/MTI 
depending on the sign of temperature gradient. The adi- 
abatic cosmic-ray instability (ACRI) described by equa- 
tion (19) is different from the MTI and HBI in that 
it cannot saturate by magnetic field reorientation (the 
cosmic-ray driving in eq. [19] does not depend on the 
field configuration term J that shows up in the ther- 
mal driving). It must lead to vigorous convection that 
changes the cosmic ray entropy profile to become nearly 
adiabatic. 

2.2. The Cosmic Ray Diffusion Coefficient 

A cosmic ray particle streaming through a magnetized 
plasma is efficiently scattered in pitch angle by magnetic 
fluctuations with wavelengths comparable to its Larmor 
radius. An unavoidable source of magnetic fluctuations 
is the self-excited streaming instability (e.g., Kulsrud & 
Pearce 1969). The effect of pitch-angle scattering due to 
Alfven waves is that the bulk speed of cosmic rays rela- 
tive to the thermal plasma is close to the Alfven speed, 
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Fig. 1. — Different timescalcs for the initial cluster model: 
isothermal sound crossing time (solid line, t sn( j = r/c s , where 
c s = [p/p] 1 / 2 ), cooling time (dotted line; we do not include cooling 
in our calculation but the cooling time is shown for comparison), 
cosmic ray energy injection time (tj n j = p/(-f — 1)Q C ; short-dashed 
line), HBI growth time (*hbi = [gdlnT/dr]~ long-dashed line), 
conduction time (short dot-dashed line; < con( j = r 2 nks / ), and 
cosmic ray diffusion time (long dot-dashed line; t^g = r 2 /D^) for 
D,| = 10 29 cnA" 1 . 

i.e., Vd — va = c 2 /uL cr , where Vd is the cosmic ray drift 
velocity relative to the thermal plasma, v is pitch-angle 
scattering rate, and L cr is cosmic ray gradient scale (Kul- 
srud 2005). The cosmic rays stream along the magnetic 
field direction, and down the cosmic ray pressure gra- 
dient. In addition to streaming with Alfven wave pack- 
ets, cosmic rays also undergo momentum-space diffusion, 
which leads to spatial diffusion along the field lines with 
Z?|| = c 2 jv = (vd — va)L ci . Thus, if cosmic ray scatter- 
ing is efficient and Vd ~ va, the diffusion timescale over 
scales comparable to L cr is much longer than the Alfven 
crossing time. 

With the above model for cosmic ray scattering one can 
self-consistently calculate the cosmic ray diffusion coef- 
ficient in terms of the plasma parameters (e.g., Loewen- 
stein, Zweibel, & Begelman 1991). In the present cal- 
culations we do not explicitly include the effects of cos- 
mic rays streaming with respect to the plasma. Instead, 
the diffusion coefficient in equation (7) should be inter- 
preted as an effective diffusion coefficient, taking into 
account both microscopic diffusion and streaming along 
turbulent magnetic field lines. The cosmic ray "diffu- 
sion" time due to cosmic rays streaming along random 
magnetic field lines can be crudely bounded by the Alfven 
crossing time (< t/va)- Together with the fact that the 
true microscopic diffusion time is much longer than the 
Alfven crossing time (as argued above), this motivates 
our choice of D\\ — arvA for the cosmic-ray diffusion co- 
efficient in most of our numerical simulations, where a 
is a factor of order unity. Somewhat arbitrarily, we take 
a = 0.4, but our results are insensitive to a so long as 
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a < 1. 

This estimate of D|| (< tva) is consistent with the mea- 
sured Galactic cosmic ray diffusion coefficient for ~ GeV 
particles, ~ 10 28 cm 2 s _1 (Berezinskii et al. 1990), using 
typical values for magnetic field strength and cosmic ray 
scale height in the Galaxy. However, at higher energies 
the diffusion coefficient increases as e ' 5 , where e is the 
cosmic ray energy (e.g., see Engelmann et al. 1990). For 
a cosmic ray energy distribution function steeper than 
e~ 2 (in the Milky Way it scales as e~ 2 - 7 from 1 to 10 5 
GeV), the cosmic ray pressure in equation (2) will be 
dominated by the lowest energy cosmic rays, and a diffu- 
sion coefficient Du < tva seems appropriate for the fluid 
description of cosmic rays considered here. With this 
choice, the ratio of the buoyancy timescale (t 2 uoy ~ r/g) 

to the cosmic ray diffusion timescale (tdis ~ r 2 /D^) is 

ibuoy/idiff < p- 1/2 . In clusters 10 < (3 < 1000 (Gov- 
oni & Ferctti 2004) so that we expect cosmic rays to be 
adiabatic on relatively large length scales and thus to be 
susceptible to cosmic ray driven convection when their 
entropy gradient is sufficiently large. Our linear stability 
analysis differs from that of Ghandran & Dennis (2006) 
and Dennis & Chandran (2008), in that we allow for a 
background heat flux. In addition, while they obtained 
the dispersion relation allowing the thermal plasma and 
cosmic rays to be either both in the diffusive limit or both 
in the adiabatic limit, we have argued that the more rel- 
evant case in clusters is likely that of diffusive thermal 
plasma and adiabatic cosmic rays. Because of uncertain- 
ties in Du, we have carried out numerical simulations for 
different values of Du. We find that for Du = 10 28 cm 2 
s -1 , the cosmic rays are effectively adiabatic at all scales, 
and even for D\\ as high as 10 29 cm 2 s -1 , the cosmic rays 
are adiabatic on large scales, r > 10 kpc (see sec. 3.2.4). 

3. NUMERICAL SIMULATIONS 

We have extended the ZEUS-MP MHD code (Hayes et 
al. 2006; Stone & Norman 1992a, b) to include thermal 
conduction along magnetic field lines (Sharma & Ham- 
mctt 2007), and have added cosmic rays as an additional 
fluid diffusing along magnetic field lines. We numerically 
solve equations (l)-(7). As mentioned earlier, we do not 
include plasma cooling. The cosmic ray energy equa- 
tion, thermal conduction in equation (4), and the cosmic 
ray pressure gradient in the equation of motion are im- 
plemented in an operator split fashion with appropriate 
source and transport terms. Both thermal conduction 
and cosmic ray diffusion are sub-cycled. 

We have tested the code extensively. The anisotropic 
cosmic ray diffusion equation is analogous to anisotropic 
thermal conduction. We use the method of Sharma & 
Hammett (2007) which preserves positivity of p cv . Ap- 
pendix shows a 1-D shock tube test, adapted from Pfrom- 
mer et al. (2006), which shows that adiabatic evolution 
of cosmic rays is accurate. The thermal conductivity is 
chosen to be the Spitzer value, 

= L84x10^ t5/2 -i K -7/2 cm -i (2Q) 
11 In A 

Based on the discussion in §2.2 we choose Du = OArvA 
for most of our calculations, where va is the local Alfven 
speed; to test the dependence on , we also carry out 



simulations with a constant value of Dm = 10 28 , 10 29 

2—1 

cm s . 

How cosmic rays are produced and distributed in the 
ICM is still pooriy understood (e.g., Eilck 2003). Thus, 
we use a simple phenomenological source term to drive 
the cosmic ray pressure in the inner parts of the cluster. 
The cosmic ray energy source term in equation (5) is 
based on Guo & Oh (2008), 

«.~^(£P ['-'-""•>']■ <»> 

We take ro = 20 kpc and v = 1.5, which leads to a cen- 
trally peaked cosmic ray entropy, and take the cosmic 
ray energy injection rate to be J Q C A-Kr 2 dr = 5 x 10 42 erg 
s^ 1 . 3 Physically, the cosmic ray energy injection rate is 
presumably related to the accretion rate onto the central 
black hole, via J Q c 4irr 2 dr = eMc 2 , where e is the effi- 
ciency of cosmic ray energy production. Indeed, Allen et 
al. (2006) show that the mechanical luminosity of jets 
can be ~ few % of the inferred Bondi luminosity. For 
M = 0.1 Mq yr _1 , our cosmic ray energy injection rate 
corresponds to £ ~ 10~ 3 and thus the level of cosmic- 
ray power used here is observationally and theoretically 
quite reasonable. 

The cosmic-ray pressure fraction and its dependence on 
radius are not that well-constrained observationally. For 
a few clusters (e.g., Virgo, Perseus, and Fornax) obser- 
vational constraints indicate that p C r/p ^ 0.2 averaged 
over the cluster core (Pfrommer & Enfilin 2004; Chu- 
razov et al. 2008). It is possible, however, that larger 
cosmic-ray pressure fractions arise in clusters in which 
AGN feedback is expected to play a particularly strong 
role, such as Hydra A or Scrsic 159-03 (see, e.g., Za- 
kamska & Narayan 2003; Chandran & Rasera 2007); 
cosmic-rays may also be particularly important at small 
radii, and quite likely contribute significantly to the pres- 
sure in and around radio bubbles (which are associated 
with a deficit of X-ray emission from the thermal plasma; 
e.g., Birzan et al. 2004). In our simulations, we choose 
parameters in equation (21) such that the cosmic ray 
pressure is smaller than the plasma pressure even at late 
times. For a larger Q c , we find that the cosmic ray pres- 
sure builds up faster than the rate at which cosmic ray 
pressure can be transported outwards via convection (or 
diffusion, but the latter is slow), and the cosmic ray pres- 
sure can become larger than the plasma pressure. In cal- 
culations with cooling, it is likely that a larger cosmic ray 
injection rate could remain consistent with p C r/p < 0.25: 
if the gas is allowed to cool then a large part of the cos- 
mic ray energy may be channeled into plasma heating 
(and thus cooling) without building up a larger cosmic 
ray pressure. 

To assess how effectively turbulence generated by the 
HBI or ACRI mixes the plasma, we solve for the advec- 
tion of a passive scalar density (e.g., a proxy for metal- 

3 A physically more realistic model would be to include a feed- 
back source term for cosmic rays where, instead of a fixed cosmic 
ray luminosity, a fixed fraction of the instantaneous mass accre- 
tion rate is converted into cosmic ray power; this level of detail 
is unnecessary for studying the basic physics of cosmic-ray driven 
convection but will be included in future calculations with radiative 
cooling. 
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licity) /, using 



(22) 



Appendix A shows the behavior of the passive scalar den- 
sity in a shock tube test. 

3.1. Simulation Parameters 

The simulations are carried out in spherical (r,0,<j>) ge- 
ometry with the inner boundary at r- m = 1 kpc and the 
outer boundary at r out = 200 kpc. Strict outflow bound- 
ary conditions are applied to the radial velocity at the 
inner and outer radial boundaries. The plasma pressure 
and density are held fixed at the outer boundary to pre- 
vent spurious oscillations; the plasma cooling time at r ou t 
is longer than the Hubble time. All other plasma and cos- 
mic ray variables are copied on ghost zones at both the 
inner and outer radial boundaries. Reflective boundary 
conditions are applied at the 9 boundaries (9 — 0, n). 
Periodic boundary conditions are applied in the (f> direc- 
tion. A logarithmic grid is chosen in the radial direc- 
tion, while the grid is uniform in 9 and 4>. Our fidu- 
cial run uses a 128 x 64 x 32 grid, with < 9 < n 
and < 4> < 2-k. With these choices, Ar/r = 0.042, 
A6 = 0.05, and A</> = 0.2. We also carried out a num- 
ber of 128 x 64 2-D (axisymmetric) simulations, and one 
higher resolution, 256 x 128, 2-D simulation for conver- 
gence studies. 

We typically initialize a weak (/3 > 10 6 everywhere) 
split-monopole magnetic field with B oc r~ 2 , although 
in one calculation (CRM), we use a monopole field to 
compare simulations with and without net magnetic flux. 
In runs with cosmic rays, the initial cosmic ray pressure 
is small (0.005 times the plasma pressure at n n ) and 
varies as r -3 ; the cosmic ray pressure builds up in time 
via the source term in eq. [5]). The initial pressure is 
chosen such that the plasma is in dynamical equilibrium 
(d[p + p CI ] / dr = —pg). Initial thermal equilibrium is not 
imposed. 

As in Guo & Oh (2008), we use cluster parameters 
relevant for Abell 2199. The gravitational potential ($) 
is the sum of the dark matter potential (&dm), 



DM 



2GM a ln(l + r/r s 
r s r/r s 



(23) 



where Mq = 3.8 x 10 14 M Q is the characteristic dark mat- 
ter mass, r s = 390 kpc is the scale radius (Navarro, 
Frenk, & White 1997), and the potential due to the 
central cD galaxy (& c d), 



$ cD = -4irGp r 



2 \n{r/r g + ^l + (r/r g f) 



r/r g 



(24) 



where r g = 2.83 kpc, p a = 5.63 x 10 



-23 g cm 3 ^ gee 



Kelson et al. 2002 for cD galaxy NGC 6166). 

The ideal gas law p = nksT is used with np — n e p e = 
p/rrip, where p (p e ) is the mean molecular weight per 
thermal particle (electron) and n (n e ) is the total (elec- 
tron) number density. We assume a fully ionized plasma 
with hydrogen mass fraction X = 0.7, helium mass frac- 
tion Y = 0.28, such that p = 0.62 and p e — 1.18. Since 
we do not include cooling, the metallicity appears only in 





Fig. 2. — Angle averaged plasma temperature (left) and electron 
number density (right) as a function of radius at different times for 
CR (upper panel) and NCR (bottom panel). Solid lines are at 1/3, 
1, 3, and 9 Gyr; dotted line is the initial profile. The temperature 
at ~ 100 kpc increases at late times while density at ~ few kpc 
decreases with time. 



TABLE 1 

Parameters for different runs 



Label 


Dim. 


initial B 


D ll 


CR angle* 


CR* 


3-D 


split-M 


OArvA 


90 u 


NCR 


3-D 


split-M 








CRM 


3-D 


monopole* 


OAtva 


90° 


CR2D 


2-D 


split-M 


OArvA 


90° 


CR2D-dbl A 


2-D 


split-M 


OAtva 


90° 


CR28 


2-D 


split-M 


10 28 


90° 


CR29 


2-D 


split-M 


10 29 


90° 


CR30 


3-D 


split-M 


OAtva 


30° 


CR30-adt 


3-D 


split-M 


OArvA 


30° 



Half-angle around the polar axis over which CR source is applied. 
*The fiducial run. 

Although monopolar, V ■ B = everywhere, including the 
boundaries, since the origin is excluded from the computational 
domain. 

A Resolution for CR2D-dbl is 256 x 128, double that of CR2D. 
Plasma is adiabatic (km = 0) for this run. 

the conversion of pressure into plasma temperature used 
to calculate the thermal conductivity in equation (20). 

To roughly match the observations (Johnstone et al. 
2002), the initial temperature increases linearly from 1.5 
keV at r in to 4.6 keV at r out (see Fig. 2); the electron 
number density is fixed to 0.0015 cm~ 3 and the tem- 
perature is fixed to 4.6 keV at r out at all times. The 
initial density, obtained from imposing hydrostatic equi- 
librium, is quite a bit larger than the observed density 
at r; n . This is because the density obtained from hydro- 
static equilibrium is extremely sensitive to the form of 
the temperature profile, for which we use a simple linear 
fit. However, since we do not include cooling, the large 
inner density does not significantly affect our results. 

3.2. Results 

Table 1 summarizes the properties of our simulations. 
Although we list a number of calculations in Table 1, 
we focus on two 3-D simulations: CR (the fiducial run) 
and NCR. Cosmic rays are not included in NCR; NCR 
thus serves as a control run that allows us to isolate the 
effects of cosmic rays. The aim of the rest of the sim- 
ulations is to understand certain aspects of the physics 
in more detail. As described earlier, for all runs except 
CRM we initialize a split-monopole magnetic field. For 
convergence studies, we carried out a two dimensional 
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(axisymmctric) version of CR, CR2D, and compared it 
with a run with double the resolution, CR2D-dbl. All 
runs, except CR28 and CR29, use D» = OAtva] a fixed 
value for D\\ is chosen for CR28 and CR29 to test the 
influence of D\\ on our results. In runs CR30 and CR30- 
ad, the cosmic ray source term is only applied within 
30° of the pole to study angular diffusion of cosmic rays 
with (CR30) and without (CR30-ad) thermal conduction 
along field lines; in the rest of the calculations, cosmic 
ray injection is spherically symmetric. 

Our initial profiles are in dynamical equilibrium, but 
not in thermal equilibrium. The plasma will remain 
static if thermal conduction and cosmic rays are not in- 
cluded, since the magnetic field is very weak and the 
plasma is stably stratified according to the Schwarzschild 
criterion. However, in the presence of thermal conduc- 
tion, the background temperature and density change in 
time; in addition, the magnetic field lines are reoriented 
by the HBI. When the cosmic ray source term is applied, 
the cosmic ray entropy gradient can drive convection and 
mixing due to the ACRI. In this section we discuss the 
influence of these effects on the structure of the ICM. 

Figure 1 shows different timescales in the initial state: 
the isothermal sound crossing time (solid line) is the 
shortest timescale; the plasma cooling time (dotted line) 
is included for comparison with the other timescales, al- 
though cooling is not included in the simulations; the 
growth-time of the HBI (long-dashed line) varies slowly 
as a function of radius; the cosmic ray injection timescale 
(short-dashed line), i.e., the timescale for the cosmic- 
ray source term to increase the cosmic ray pressure by 
an amount comparable to the plasma pressure, increases 
rapidly with radius since the cosmic ray source term is 
centrally concentrated; the thermal conduction timescale 
(short dot-dashed line) has a maximum at intermediate 
radii and is comparable to the HBI timescale at both 
rj n and r out ; finally, the cosmic ray diffusion timescale 
(long dot-dashed line; shown for D» = 10 29 cm 2 s _1 ) is 
shorter than the buoyancy timescale only within ~ 10 
kpc. Since the diffusion time is oc 1/Dn, this implies 
that cosmic rays are effectively adiabatic for smaller D\\ , 
i.e., in all runs except CR29 (Table 1). 

3.2.1. The Fiducial Run (CR) 

We discuss the fiducial simulation (CR) in detail and 
compare it with the simulation without cosmic rays 
(NCR). Both simulations CR and NCR show similar 
properties for radii > 30 kpc, the radius outside of which 
the cosmic ray pressure is always small. Since there is 
no cooling to balance heating by thermal conduction in 
the initial state, the initial thermal properties will be 
modified on the conduction timescale. 

Figure 2 shows the angle averaged temperature and 
density profiles as a function of radius for different times 
for runs CR and NCR. Figure 2 shows that for run NCR, 
the temperature becomes isothermal near the inner and 
outer radii due to thermal conduction along initially ra- 
dial magnetic field lines. Although Figure 1 shows that 
the HBI timescale is always shorter than the conduction 
time, the two are of the same order at both ri n and r ou t 
where thermal conduction makes the plasma isothermal 
before the HBI can grow significantly. The HBI, which 
is active within 30 kpc < r < 100 kpc, reorients the 
magnetic field lines to be primarily perpendicular to the 
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Fig. 3. — Time (from tij/2 to 3i///4) and angle averaged an- 
gle between magnetic field unit vector and the radial direction 
(cos _1 (|6 ■ r\) in degrees) as a function of radius: for CR (solid 
line), for NCR (dot-dashed line), for CR2D (dotted line), and for 
CR2D-dbl (dashed line). 

radial direction (see Fig. 3). This creates a thermal 
barrier and a large temperature gradient at these radii. 
The formation of such a thermal barrier is not forced by 
the boundary condition since the temperature at r; n is 
floating; if the conduction timescale were shorter than 
the HBI timescale at all radii, the plasma would become 
isothermal (4.6 kcV) at all radii. Run CR shows similar 
behavior at larger radii but differs substantially at small 
radii. Since the cosmic ray pressure provides a substan- 
tial fraction of the pressure support to balance gravity 
within 30 kpc, the plasma density and pressure at small 
radii in run CR is substantially smaller than in NCR. 
The plasma temperature is similar in magnitude for CR 
and NCR but shows a minimum at <~ 20 kpc for CR at 
late times. 

Figure 4 shows the cosmic ray entropy profile (top) and 
the ratio of cosmic ray pressure to plasma pressure (p C r/p, 
bottom) as a function of radius for different times; we de- 
fine cosmic ray entropy as p„/p 7 " since this is the quan- 
tity whose gradient determines the buoyant response of 
adiabatic cosmic rays (see, e.g., eq. [19]). Since the cos- 
mic ray source term is chosen to be a strong function of 
radius (Q c oc r -2 - 5 for r < 20 kpc), it drives convection 
due to the ACRI when the cosmic ray pressure builds 
up and becomes comparable to the plasma pressure. At 
later times, cosmic ray injection does not continuously in- 
crease the inner cosmic ray pressure with time; instead, a 
cosmic ray driven convection front spreads radially out- 
wards. Convection drives the cosmic rays to be adia- 
batic (pcr//o 7cr ~ constant) in regions where the cosmic 
ray pressure is not negligible compared to the plasma 
pressure. 

Figure 5 shows 2-D contour plots of turbulent veloc- 
ity (absolute value of velocity) in the 4> — n plane at 9 
Gyr for CR and NCR. Velocities are similar for r > 30 




Fig. 4. — Angle averaged cosmic ray entropy (arbitrary units; 
top), and the ratio of the cosmic ray pressure to the thermal plasma 
pressure (bottom) as a function of radius for 1/3, 1, 3, and 9 Gyr 
(solid lines); dotted line is the initial profile. The ACRI flattens 
the cosmic ray entropy profile and both p C r/p 7cr and p C i/p increase 
and move outwards in time. 
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Fig. 5. — Contour plot (in the <j> = n plane) of Log\QV (in cm 
s" 1 ) at 9 Gyr for CR (left) and NCR (right). Turbulent velocities 
> 100 km s — 1 are attained in inner 20 kpc for cosmic-ray driven 
convection (left). The turbulent velocities are significantly smaller 
in the absence of cosmic rays. Arrows show the r — 9 projection of 
the velocity unit vector. 

kpc where cosmic ray pressure is negligible and fluid mo- 
tions are driven by the HBI. The maximum turbulent 
velocity in the HBI-dominated regions is ~ 30 km s . 
Run CR shows much larger turbulent velocities (<~ 100 
km s _1 ) for r < 30 kpc; note that the turbulent veloc- 
ities induced by the ACRI are consistent with mixing 
length theory, with L c ~ Airr 2 pv\, where L c is the power 
carried by convection and v c is the resulting convective 
velocity. The large turbulent velocities in the presence of 



the ACRI are sufficient to prevent catastrophic cooling 
according to the models of Chandran & Rasera (2007). 
By contrast, the turbulent velocities are extremely small 
at r < 30 kpc for NCR because the plasma temperature 
gradient is wiped out by conduction before the HBI can 
drive any turbulence. 

While the turbulent velocity vectors are roughly 
isotropic at r < 30 kpc for CR, they are aligned primar- 
ily perpendicular to the radial direction at large radii 
where the HBI dominates. This is consistent with the 
result that while the HBI saturates by reorienting mag- 
netic field lines perpendicular to gravity (e.g., Parrish 
& Quataert 2008), the ACRI drives roughly isotropic 
convection irrespective of the magnetic field geometry. 

Figure 3 shows the time-averaged (from i# /2 to 3t jj/4) 
angle between the magnetic field vector and the radial 
direction as a function of radius for runs CR and NCR; 
the angle is defined with respect to the radial direction, 
such that it is 0° for the initial split-monopole field. The 
average angle is similar for r > 30 kpc, where the cos- 
mic ray pressure is negligible, for both simulations. For 
run CR, radii r < 30 kpc are convectively stirred by 
the ACRI and the average angle between the magnetic 
field and the radial direction is w 55°, close to the value 
expected for a uniform, random magnetic field unit vec- 
tor (cos- 1 1/2 = 60°). For both CR and NCR at 30 kpc 
< r < 100 kpc, the magnetic field is nearly perpendicular 
to the radial direction because of the HBI. The average 
angle between the magnetic field unit vector and the ra- 
dial direction at these radii is w 75°. For r > 100 kpc 
and for r < 20 kpc in run NCR, the HBI is weak since 
thermal conduction wipes out the temperature gradient 
before the HBI can grow significantly; as a result, the 
field lines are not perpendicular to the radial direction. 

Figure 6 shows angle and time averaged (from ijj/2 to 
3iff/4) kinetic and magnetic energy profiles as a func- 
tion of radius for runs CR, CRM, and NCR. Both the 
magnetic and kinetic energies are generally amplified, al- 
though to varying degrees. The kinetic energy is very 
effectively amplified at small radii by the action of the 
ACRI in run CR; the turbulent Mach number is ~ 0.1 
(see also Fig. 5). The magnetic energy as a function of 
radius in CR is quite striking in that there is no magnetic 
energy enhancement at small radii: there is a slight bump 
in magnetic energy at 20 kpc but the final magnetic en- 
ergy is smaller than the initial magnetic energy for r < 3 
kpc. In turbulent dynamos one often finds the turbulent 
magnetic energy to be of the order of turbulent kinetic 
energy (e.g., Cho & Vishniac 2000); this is clearly not 
the case for the run CR. Magnetic field amplification can, 
however, be subtle; e.g., in local shearing box simulations 
of the magnetorotational instability with no net magnetic 
flux, magnetic field amplification (and associated stress 
and turbulence) occurs only for Prandtl numbers exceed- 
ing unity (e.g., Fromang et al. 2007). We do not have 
explicit viscosity and resistivity and it is possible that 
the effective Prandtl number in our simulation is small, 
so that dissipation of the initially split monopolar field 
at the grid scale dominates over magnetic field enhance- 
ment by convective turbulence. To better understand 
this, we have done a simulation with cosmic rays with a 
net initial magnetic flux (run CRM), but with all other 
properties of the simulation the same. Figure 6 shows 
that in this case, the magnetic energy density is much 
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r (kpc) 

Fig. 6. — Time (from to 34j^/4) and angle averaged mag- 

netic and kinetic energy densities (in erg cm -3 ) as a function of 
radius for CR (solid line), CRM (long dashed line), and NCR (dot- 
ted line). The kinetic energy density is larger than the magnetic 
energy density in all cases. The initial magnetic energy density pro- 
file (short dashed line) is also shown for comparison. Run CRM, 
which includes a mean field, results in a much larger amplification 
of the magnetic energy as compared to CR. 

larger in the inner regions as compared to CR, although 
the magnetic energy is still ~ 100 times smaller than the 
kinetic energy. In this case, the magnetic field strength 
in the center of the cluster is amplified to ~ 0.1-1 /iG by 
convective motions. The kinetic energy density profiles 
are almost identical for CR and CRM, indicating that 
the properties of the convection are not very different in 
the two cases, although the efficiency of magnetic field 
amplification differs dramatically. 

At larger radii (30 kpc < r < 100 kpc), the HBI causes 
amplification of the magnetic and kinetic energies in both 
CR and NCR. The magnetic energy is enhanced by a fac- 
tor ~ 100 at r ss 60 kpc. This level of field amplification 

- a factor of ~ 10, primarily of the 8 and <j) components 

- is what is required to reorient initially radial magnetic 
fields into fields that are perpendicular to the radial di- 
rection (as was seen in previous HBI and MTI simula- 
tions; Parrish & Quataert 2008; Parrish & Stone 2007; 
Sharma, Quataert, & Stone 2008). The kinetic energy is 
also enhanced at these radii because of turbulence driven 
by the HBI. 

In order to study the numerical convergence of our fidu- 
cial run, we carried out an axisymmetric 2-D run anal- 
ogous to CR - CR2D - and an axisymmetric run with 
double the resolution - CR2D-dbl (see Table 1); study- 
ing convergence directly with the 3D run CR would have 
been computationally prohibitive, requiring ~ 32 times 
more cpu time. The results from runs CR2D and CR2D- 
dbl are nearly identical to each other; in particular, an- 
gle averaged plots such as those shown in Figures 2, 4, 
3, and 6 are quite similar in the two cases. We have not 
included these 2-D results in every figure, but to illus- 
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Fig. 7. — Contour plot (at <f> = it) of the passive scalar den- 
sity, / (normalized to the initial maximum), at 9 Gyr for CR (left) 
and NCR (right). The passive scalar is initialized near the origin. 
Projection of magnetic field unit vector is over-plotted by arrows. 
While the passive scalar diffusion is negligible for NCR, turbulent 
mixing is efficient within 30 kpc for the cosmic-ray driven convec- 
tion (left). 

trate the basic convergence result, Figure 3 shows that 
the average angle between the magnetic field unit vector 
and the radial direction is quite similar for CR2D and 
CR2D-dbl (which are both similar to the 3-D run CR). 
It is interesting to note that the 2-D runs differ from the 
3-D run CR in one important respect: the amplification 
of the magnetic field at small radii is significantly larger 
in 2-D than in 3-D; the magnetic field energy density 
at small radii in the 2-D runs is comparable to the 3-D 
run with a net magnetic flux (CRM) in Figure 6. 4 As 
mentioned earlier, the reason for the lack of significant 
magnetic field amplification in the run CR is not entirely 
clear. 

3.2.2. Diffusion of a Passive Scalar 

Figure 7 shows /, the passive scalar density, in the in- 
ner 40 kpc for CR and NCR at 9 Gyr, as well as the 
projection of the magnetic field unit vectors. The pas- 
sive scalar density is initialized to be a large number 
(/ = 10 15 ) for r < 1.25 kpc (corresponding to four radial 
zones) and is negligible (/ = 10~ 15 ) for r > 1.25 kpc. 
The goal of initializing a passive scalar is to study mix- 
ing due to turbulence. Observations of clusters reveal a 
metallicity distribution that is more spatially extended 
than the light distribution of the central galaxy. This 
may indicate turbulent transport of metals in clusters 
(e.g., Rebusco et al. 2006; Rasera & Chandran 2008). 
As expected, run CR with large turbulent velocities at 
r < 30 kpc also results in efficient mixing. Mixing is 
negligible for NCR because the inner radii (r < 30 kpc) 
are isothermal and are thus not stirred by the HBI. For 
a direct comparison with observations, one must include 
a time dependent, spatially distributed source term in 
the passive scalar equation which represents metal en- 
richment due to Type la supernovae (e.g., Rebusco et 

4 By the anti-dynamo theorem, the amplification in 2D must be 
transient. The dynamical time in clusters is so long, however, that 
the "transient" can last a Hubble time! 
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Fig. 8. — Angle averaged passive scalar density (Logiof ; normal- 
ized to the initial maximum) as a function of radius for CR at 1/3, 
1, 3, 9 Gyr (solid lines); the initial profile, shown by the dotted 
line, very close to the y axis. Passive scalar density decreases as 
it spreads out with time. For comparison, the dashed line shows 
a Gaussian fit with a diffusion coefficient D = 10 28 cm 2 s — 1 at 9 
Gyr (~ exp[-r 2 /6Dt]). 

al. 2005); this is beyond the scope of the present paper. 
Nonetheless, our results indicate that cosmic-ray driven 
convection is an efficient mechanism for mixing plasma 
in clusters of galaxies. 

Figure 8 shows the angle averaged passive scalar den- 
sity as a function of radius for run CR at 1/3, 1, 3, 9 Gyr. 
Also shown is a Gaussian fit (at 9 Gyr) with a diffusion 
coefficient of 10 28 cm 2 s _1 . The passive scalar density 
at 9 Gyr is flatter than the Gaussian fit at < 20 kpc, 
implying that diffusion due to convection driven by cos- 
mic rays corresponds to an effective diffusion coefficient 
somewhat larger than 10 28 cm 2 s _1 . Beyond ~ 30 kpc, 
the cosmic ray pressure is unimportant, and turbulence 
is driven by the HBI. This change in the source of the 
turbulence accounts for the rapid decline in the passive 
scalar density at large radii in Figure 8. 

To isolate just the mixing induced by turbulence driven 
by the HBI, we have taken run NCR at t H /4(= 3.425 
Gyr) and initialized a passive scalar density peaked at 
r w 53 kpc, the radius where the temperature gradient 
is large and the HBI is active (see Fig. 2). The pas- 
sive scalar is initialized to be 10 15 at two grid points at 
r = 53 kpc, 9 = 7r/2 and (j) = tt, and negligible (10~ 15 ) 
elsewhere. Figure 9 shows if> = tt snapshots of the passive 
scalar density at later times. The passive scalar diffuses 
more rapidly in the 9 direction. This is because the tur- 
bulent velocities due to the HBI are larger in the direc- 
tion perpendicular to gravity than they are in the radial 
direction (just as the magnetic field components perpen- 
dicular to gravity are preferentially amplified). To esti- 
mate the diffusion coefficient in the r and 9 directions, we 
compare how much the passive scalar has spread in the 
two directions; the full width at half maximum (FWHM) 
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Fig. 9. — Contour plot (at <f> = 7r) of the passive scalar density 
(/, normalized to the initial maximum) at 3.7, 4.25, 5.9, 10.8 Gyr 
for the run NCR. The passive scalar density is initialized in a small 
region (two grid points along r, 8, <j>) centered at 53 kpc at 3.425 
Gyr, and is negligible everywhere else. Turbulent mixing due to the 
HBI is faster in the 8 direction compared to the radial direction. 

for / along r and 9 at 9 Gyr is w 10, 40 kpc, respectively. 
For comparison, the FWHM at 9 Gyr for / for run CR 
shown in Fig. 8 is ~ 60 kpc. Thus the diffusion co- 
efficient due to the HBI alone is ~ 2 (perpendicular to 
gravity) and ~ 50 (parallel to gravity) times smaller than 
the diffusion coefficient due to the ACRI in run CR. Al- 
though these precise numerical values likely depend on 
the detailed parameters of our simulations, the fact that 
the HBI primarily induces turbulence and mixing in the 
plane perpendicular to gravity is a generic result. 

3.2.3. Heat Flux modified by the HBI and ACRI 

Many 1-D models of clusters parameterize thermal con- 
duction by its ratio to the Spitzer value (Zakamska & 
Narayan 2003; Chandran & Rasera 2007; Guo & Oh 
2008). However, our simulations show that, because of 
plasma instabilities that operate in clusters (e.g., HBI 
and ACRI), a reduction of the conductivity by a fixed 
factor is not applicable (see also Parrish & Stone 2007; 
Parrish & Quataert 2008; Sharma, Quataert, & Stone 
2008). At large radii in cluster cores, where the cos- 
mic ray pressure is negligible, the HBI can orient field 
lines perpendicular to the radial direction, but at small 
radii where cosmic rays can be significant, the magnetic 
field may be significantly more radial. For example, for 
run CR the average angle of the magnetic field relative 
to the radial direction for r < 30 kpc is ~ 55° (sec 
Fig. 3), corresponding (roughly) to a reduction factor 
of (6 • f) 2 w 1/3. For 30 kpc < r < 100 kpc, however, 
the turbulence is dominated by the HBI, and the average 
angle between the magnetic field and the radial direction 
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Fig. 10. — Angle averaged cosmic ray entropy (upper panel; arbi- 
trary normalization) and the ratio of cosmic ray to plasma pressure 
(lower panel) as a function of radius for the runs CR28 (left) and 
CR29 (right), which use a fixed cosmic ray diffusion coefficients of 
On = 10 28 and 10 29 cm 2 s _1 , respectively. The solid lines show 
profiles at 1/3, 1, 3, 9 Gyr, with the entropy and pressure ratio 
increasing in time; the initial profiles are shown with dotted lines. 
Profiles for CR28 look similar to profiles for CR in Figure 4. The 
profiles for CR29 are different in that p C r/p 7cr and p C r/p increase 
towards a maximum at the intermediate radii; this is because the 
cosmic-rays are no longer adiabatic for large Dy . 

is « 75°, corresponding to a reduction factor of « 0.07. 

An even more subtle result is that the HBI can change 
the background temperature gradient by forming ther- 
mal barriers (see Fig. 2), which will be absent with 
isotropic conduction. For example, the temperature gra- 
dient at 30 kpc < r < 100 kpc for CR and NCR at late 
times is 3 times larger than the initial temperature 
gradient. Thus, the HBI not only reduces the conductive 
heating by a factor of (6 • f) 2 « 0.07, it also increases it 
by making the temperature gradient larger by a factor 
w 3. Approximating the conductivity of a magnetized 
plasma by a constant factor with respect to the Spitzcr 
value misses all of this interesting dynamics. Whether 
this is important in real clusters with radiative cooling 
and various sources of heating remains to be seen. 

3.2.4. Runs with larger Dy (CR28 & CR29) 

We have also carried out 2-D simulations with larger 
cosmic ray diffusion coefficients since the value of 

the cosmic ray diffusion coefficient is poorly constrained 
(see §2.2). Run CR28 uses D\\ = 10 28 cm 2 s"\ the 
cosmic ray diffusion coefficient estimated for GeV cosmic 
rays in the Galaxy. Run CR29 uses D\\ = 10 29 cm 2 
s _1 . All other parameters and initial conditions are same 
as the fiducial run. Since we are comparing these 2-D 
simulations with CR, which is a 3-D simulation, we have 
verified that the run CR2D gives results similar to the 
3-D results presented here; in particular, the profiles for 
Per/ P 1 " and Pcr/p are identical to the profiles for CR 
shown in Figure 4. 

Figure 10 shows the angle averaged cosmic ray entropy 
(Pcr/p* 1 ") and the ratio of the cosmic ray pressure to 
the plasma pressure (p C r/p) as a function of radius for 
CR28 and CR29. The profiles for CR28 and CR (see 
Fig. 4) look similar; cosmic ray entropy and pressure 
are slightly more radially spread out for CR28. The 
profiles for CR29 are quite different. The entropy and 
pressure ratio profiles have a peak at ~ 20 kpc; this is 
the radius beyond which the cosmic ray diffusion time is 
longer than the buoyancy timescale (see Fig. 1). The 
ratio Pcr/p is smaller in CR29 as cosmic rays are spread 
out over a larger volume by diffusion. The cosmic ray 



entropy increases outwards for r < 20 kpc since cosmic 
ray diffusion, and not convection driven by cosmic rays, 
dominates the outward cosmic ray transport. This is dif- 
ferent from CR and CR28 where cosmic ray diffusion is 
sub-dominant. Even in CR29, cosmic rays are effectively 
adiabatic for r > 20 kpc; cosmic-ray driven convection 
is absent at these radii, however, because the cosmic ray 
pressure is not large enough to drive the ACRI (see Fig. 
10). At smaller radii, the cosmic rays are nearly isobaric 
because of rapid diffusion, and the system is formally 
unstable to the CR mediated version of the MTI (cq. 
[18]). However, because the field lines are nearly radial 
at these radii, the growth-rate of the CRMTI is quite 
slow and we do not see any indications that it develops 
in our simulations. 

To explicitly study the possibility of the ACRI setting 
in at larger radii in the cluster core, we carried out a 
simulation with the cosmic ray source term (eq. [21]) 
three times larger than in run CR29. This larger source 
term increases the CR pressure at large radii, and at late 
times Pcr/p is large enough to drive the ACRI. The tur- 
bulent velocities are ~ 100 km s _1 at r ~ 20 — 30 kpc, 
where the cosmic rays arc effectively adiabatic in spite 
of the large Du. Thus, even in the presence of rapid 
cosmic ray diffusion, the ACRI can set in at large radii 
where the cosmic rays are adiabatic, provided the cosmic 
ray pressure is sufficiently large; this may naturally oc- 
cur in the vicinity of cosmic-ray filled buoyant bubbles. 
More generally, our results demonstrate that so long as 
-^11 1 — 3 x 10 28 , cm 2 s _1 , cosmic rays will behave 
effectively adiabatically throughout the cluster core and 
bulk transport by convection and other mechanisms will 
dominate the diffusive transport. 

3.2.5. Runs with Cosmic ray Sources at the Poles (CR30 & 

CR30-ad) 

It is very unlikely that cosmic rays in clusters are in- 
jected spherically symmetrically. Instead, the injection 
likely occurs preferentially in the polar direction. To 
study the resulting physics in this case, we carried out 
3-D simulations in which the cosmic ray source term is 
applied only within 30° of the pole: CR30 and CR30-ad. 
Except for this difference all parameters for run CR30 are 
the same as run CR. Run CR30-ad differs from CR30 in 
that the plasma is adiabatic, i.e., thermal conduction is 
not included. One of the aims of these simulations is to 
show the dramatic differences that result from including 
anisotropic thermal conduction (relative to a more typi- 
cal adiabatic simulation). Cluster plasmas are observed 
to be stable to adiabatic convection because the entropy 
increases outwards (e.g., Piffaretti et al. 2005). How- 
ever, convection in an anisotropically conducting plasma 
depends on the temperature gradient, and not the en- 
tropy gradient, and the system is unstable independent 
of the sign of the temperature gradient. This makes it 
much easier to mix a thermally conducting plasma than 
an adiabatic plasma. In clusters, this implies that turbu- 
lence produced by external means, e.g., the ACRI, wakes 
due to galaxy clusters, etc., may be an effective way of 
mixing the thermal plasma. 

Figure 11 shows contour plots (cf> = tt snapshot) of the 
ratio of cosmic ray pressure to plasma pressure (p C r/p) at 
1/3, 1, 3, 6 Gyr, for CR30 and CR30-ad. For CR30, the 
cosmic rays become unstable to the ACRI in the polar 
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Fig. 11. — Contour plots within 40 kpc, in the <f> = tt plane, 
of Logio(p C r/p) for runs CR30 (top) and CR30-ad (bottom) at 
t = 1/3, 1, 3, and 6 Gyr (from left to right). The ratio 
Pci/p is not shown if it is smaller than 10 — 4 . The adiabatic 
plasma in CR30-ad (bottom) artificially suppresses the angular 
and radial mixing of the relativistic and thermal plasma that 
is present in the simulations with anisotropic thermal conduc- 
tion (CR30; top). For movies corresponding to this figure, see: 
http: / /astro. berkeley.edu/ ~psharma/clustermovie.html. 

region. The resulting turbulence is able to convectively 
mix the plasma, not only in the unstable radial direction, 
but also in the marginally stable 9 direction. Instead of 
cosmic rays being confined only to the 9 — 30° cone, con- 
vection effectively mixes plasma in both the radial and 
angular directions. In addition, at radii beyond 20 kpc 
where the temperature gradient is appreciable (see Fig. 
2), the HBI mixes material primarily in the 9 direction, 
as seen by the 9— oriented fingers at late times in Figure 
11. For CR30, p C r/p ~ 0.02 at 3 Gyr in the equatorial 
region for r < 10 — 20 kpc. Although cosmic rays are 
still dominant near the pole, convection brings a non- 
negligible amount cosmic rays into the equatorial region. 
In comparison, there is no sign of convective overshoot 
in CR30-ad because convective motions in the thermal 
plasma near the equator are strongly stabilized by a large 
positive plasma entropy gradient. The polar cosmic ray 
dominated plasma does expand somewhat in both r and 
9 as it becomes over-pressured. However, the value of 
Pcr/p at 3 Gyr is < 10 everywhere in the equatorial 
plane for CR30-ad. 

4. SUMMARY & ASTROPHYSICAL IMPLICATIONS 

The X-ray emitting plasma in clusters of galaxies is 
hot (T ~ 1 - 10 keV) and dilute (n - 0.001 - 0.5 cm" 3 ), 
so that the transport of heat and momentum along mag- 
netic field lines can be energetically and dynamically im- 
portant. In addition, jets launched by a central AGN 
produce relativistic plasma (cosmic rays), which are ob- 
served in part as bubbles of radio emission associated 
with deficits of thermal X-ray emission (X-ray cavities; 
e.g., Birzan et al. 2004). In this paper we have studied 
the transport properties of an ICM composed of cosmic 
rays and thermal plasma. We have argued that cosmic 
ray diffusion is likely to be slow because of scattering by 
self-generated Alfven waves (§2.2); as a result, the cos- 



mic rays are adiabatic on moderately large length scales 
> 1-10 kpc. More concretely, cosmic rays are adiabatic 
on the scale of cluster cores, so long as their parallel dif- 
fusion coefficient satisfies Du < 10 cm 2 s -1 . 5 

It is now well established that anisotropic conduction 
and anisotropic cosmic ray diffusion can dramatically 
modify buoyancy instabilities in low collisionality sys- 
tems, producing qualitatively new instabilities such as 
the MTI, HBI, and their cosmic ray counterparts (e.g., 
Balbus 2000; Chandran & Dennis 2006; Parrish & Stone 
2007; Quataert 2008; Parrish & Quataert 2008). Non- 
linear studies of these instabilities (including those in 
this paper) have demonstrated that they saturate by ap- 
proaching a state of marginal stability to linear pertur- 
bations, just as in hydrodynamic convection. However, 
in a magnetized plasma, there is an additional degree of 
freedom that is not present in hydrodynamic convection, 
namely the local direction of the magnetic field. The 
primary mechanism by which these diffusive buoyancy 
instabilities saturate is by rearranging the magnetic field 
lines, so that the linear growth rate becomes extremely 
small (see Fig. 3). This is different from entropy gradi- 
ent driven convection in adiabatic fluids, which saturates 
by producing convection that wipes out strong entropy 
gradients. Even nonlinearly, most of the energy flux in 
systems unstable to the MTI and HBI is transported by 
thermal conduction, rather than convection. Moreover, 
the saturation of these instabilities is quasilinear in the 
sense that the saturated magnetic energy is proportional 
to the initial magnetic field energy (Sharma, Quataert, 
& Stone 2008). 

In this paper we have shown analytically and through 
numerical simulations that when cosmic rays have appre- 
ciable pressure, p C r/p ^ 0.25, and an outwardly decreas- 
ing entropy (p cr /p 7cl ), they can drive strong convection 
and mixing in clusters of galaxies. This adiabatic cosmic 
ray instability ( ACRI) in the central regions of clusters of 
galaxies is a cosmic ray analogue of hydrodynamic con- 
vection familiar in the context of stars and planets. In 
particular, the nonlinear saturation of adiabatic cosmic 
ray convection is similar to that of hydrodynamic con- 
vection, and thus quite different from the saturation of 
the MTI and HBI. Our simulations of cluster cores also 
provide insight into the global saturation of the HBI at 
radii in clusters where the cosmic ray pressure is neg- 
ligible (~ 30 — 100 kpc in our models), and thus where 
the only convective instability is that driven by the back- 
ground conductive heat flux. More specifically, the pri- 
mary results of this paper include: 

• If the cosmic ray entropy decreases outwards 
and if p C r/p ^ 0.25, convection driven by the 
ACRI sets in. In the saturated state, the cosmic 
ray entropy profile becomes nearly constant 
in the region with significant cosmic ray pres- 
sure (Fig. 4). The resulting turbulent velocities 
are consistent with mixing length theory, with v c ~ 

5 A purely thermal electron-ion plasma can also show an adi- 
abatic, rather than diffusive, response even in the presence of 
rapid electron thermal conduction; this occurs if electrons and 
protons are not collisionally coupled on the buoyancy timescale. 
We find, however, that even in the outer parts of clusters, the 
electron-proton energy exchange time is shorter than the buoyancy 
timescale and thus the MTI/HBI limits are appropriate. 
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100(L C /10 43 ergs" 1 ) 1 / 3 (n/0.1 cm- 3 )- 1 / 3 (r /20kpc)- 
km s _1 where L c is the total power supplied to 
cosmic rays and tq is the pressure height scale 
of cosmic rays. The ACRI generates turbulent 
motions more effectively in cluster cores than the 
HBI alone. 

• The ACRI drives roughly isotropic convection with 
the average angle between the held lines and the 
radial direction ~ 55°; by contrast, the HBI gener- 
ates magnetic field lines that are primarily in the 
9 and <f> directions (Fig. 3), shutting off the radial 
conduction of heat. The effective radial conductiv- 
ity of a cluster plasma thus depends sensitively on 
which of these instabilities operates at a given lo- 
cation, and may not be adequately approximated 
as a fixed fraction of the Spitzer value throughout 
the cluster. 

• We have quantified the mixing of a passive scalar 
by the ACRI and HBI: the ACRI produces roughly 
isotropic mixing with a turbulent diffusion coeffi- 
cient D > 10 28 cm 2 s _1 (Fig. 7); mixing length 

1 /3 

theory predicts that D oc v c oc Lc'. At lar ger 
radii, only the HBI operates and the mixing is pri- 
marily in the 9 and <f> directions, rather than in the 
radial direction (Fig. 9). Both the ACRI and the 
HBI may contribute to mixing metals in clusters 
by redistributing, in both radius and angle, metals 
produced by Type la supernovae. Some observa- 
tions of metallicity gradients in clusters have in- 
ferred mixing at levels comparable to those found 
here (e.g., Rebusco et al. 2006). 

• It is considerably easier to mix thermal plasma 
in the presence of anisotropic thermal conduction, 
since the plasma is formally always buoyantly un- 
stable and thus already prone to mixing! By con- 
trast, treating the plasma as adiabatic (i.e., ignor- 
ing thermal conduction) results in an artificially 
stabilizing entropy gradient in cluster plasmas. 6 As 
a concrete example of these effects, we have demon- 
strated that cosmic rays initially injected into the 
polar regions can be partially mixed to the equator 
by convective overshooting in the ACRI and HBI 
unstable regions (Fig. 11); this effect is largely ab- 
sent in simulations that treat the plasma as adia- 
batic. If wave heating due to cosmic ray streaming 
or heating due to Coulomb interactions is impor- 
tant in clusters (e.g., Guo & Oh 2008), a mecha- 
nism similar to that described here may be crucial 
in redistributing cosmic rays throughout the clus- 
ter volume. More generally, to study the mixing 
produced by external sources of turbulence such as 
galactic wakes or cosmic-ray filled bubbles, we sus- 
pect that anisotropic thermal conduction must be 
accounted for, so that the buoyant response of the 
thermal plasma is correctly represented. 

Having summarized our primary results, we now de- 
scribe several caveats and directions for future research. 



2 / 3 First, to generate the ACRI, we have injected cosmic 
rays using a subsonic source term at small radii. In re- 
ality a significant fraction of the cosmic rays produced 
by an AGN are expected to be produced in a supersonic 
jet shocking against the ICM. The spatial distribution 
of cosmic rays produced by jets is poorly understood. 
The intuition drawn from our simulations should apply 
as long the source of cosmic rays ultimately produces a 
centrally concentrated bubble of relativistic plasma that 
expands subsonically. 

Our calculations intentionally do not include plasma 
cooling. Instead of trying to solve the cooling flow prob- 
lem, our goal has been to study the basic physics of buoy- 
ancy instabilities in the combined relativistic + thermal 
plasma, implicitly assuming that some heating process 
is preventing catastrophic cooling of the plasma. Our 
calculations also do not include anisotropic ion viscos- 
ity which is w 40 times smaller than electron thermal 
conductivity. Finally, we do not treat the effects of cos- 
mic ray streaming with respect to the thermal plasma 
from first principles, although our choice of the cosmic 
ray diffusion coefficient qualitatively accounts for limits 
on cosmic ray streaming produced by self-excited Alfven 
waves (see §2.2). In future work, we intend to include all 
of the above effects, which will provide a more quantita- 
tive model of plasma in cluster cores. 

Finally, we note that in a full cosmological context, 
galaxy clusters will contain a large number of galaxies 
and other dark matter substructure. The motion of such 
bound objects through the ICM will reorient the mag- 
netic field and generate downstream turbulence. The 
interplay between this turbulence and that generated by 
the instabilities studied in this paper is worth investigat- 
ing in detail in future work. This interaction may create 
a magnetic dynamo in the ICM that is more effective 
than that produced by the HBI alone: galaxies moving 
through the ICM will comb out the magnetic field lines 
in the radial direction, while the HBI will amplify the 
field and generate a strong perpendicular magnetic field 
component from the seed radial field created by galactic 
wakes. 
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6 Even including isotropic thermal conduction reduces the stabi- 
lizing effect of the entropy gradient; it does not, however, capture 
the MTI/HBI, which are driven by anisotropic thermal conduction 



along magnetic field lines. 
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APPENDIX 
NUMERICAL TESTS 
Cosmic ray shock tube 

We have tested the adiabatic implementation of cosmic rays with a 1-D shock tube problem discussed in Pfrommer 
et al. (2006) and in Rasera & Chandran (2008). Like Rasera & Chandran (2008), we also use 1024 grid points. 
Thermal conduction and cosmic ray diffusion are absent for this test problem. The left (1 < x < 1.5) and right 
(1.5 < x < 2) states are given by (p L , vl, Pl, Pctl) = (1, 0, 6.7 x 10 4 , 1.3 x 10 5 ) and (p R , v R , p R , p crR ) = (0.2, 0, 
2.4 x 10 2 , 2.4 x 10 2 ), respectively. Figure 12 shows profiles at t = (solid line) and at t = 4.4 x 10~ 4 (shorter than the 
crossing time; points). The profiles match very well with the analytic result and with Fig. 3 of Rasera & Chandran 
(2008). In addition we also test advection of a passive scalar density governed by equation (22). The passive scalar 
density shows a discontinuity at the location of the contact discontinuity; volume integrated / is not conserved but 
volume integrated pf is. 
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Fig. 12. — Profiles for different fluid variables for the shock tube test: initial profiles (solid line) and profiles at t = 4.4 X 1CP 4 (points). 
While the shock is resolved by 4 points, contact discontinuity requires more grid points to be resolved (this is a feature of all methods that 
do not solve the full Riemann problem). 



